% This 
clear all; close all;
%% 1. Model setup: please copy and paste the model set up code section 1
ModelSettings.SiteNum       = 100;       % Number of sites
ModelSettings.LengthUnit    = 'layer';   % 
ModelSettings.SiteWidth     = 2/sqrt(3); % The round particle

ModelSettings.x_max_EW      = pi;        % The x domain in the EW equation ranges from -pi to pi
ModelSettings.x_min_EW      = -pi;
ModelSettings.dx_EW = (ModelSettings.x_max_EW-ModelSettings.x_min_EW)/(2*ModelSettings.SiteNum);
dx_EW = ModelSettings.dx_EW;
x_EW = ModelSettings.x_min_EW:dx_EW:ModelSettings.x_max_EW;

ModelSettings.x_max         = ModelSettings.SiteNum*ModelSettings.SiteWidth;
ModelSettings.x_min         = 0;
ModelSettings.dx = (ModelSettings.x_max-ModelSettings.x_min)/(2*ModelSettings.SiteNum);
dx = ModelSettings.dx;
x = ModelSettings.x_min:dx:ModelSettings.x_max;

ModelSettings.mode          = 200;

K_alpha = zeros(ModelSettings.mode,1);
K_beta  = zeros(ModelSettings.mode,1);
for n=1:ModelSettings.mode
    h_sin = sin(n*x_EW);
    dh_sin = h_sin;
    dh_sin(1) = h_sin(1)-h_sin(end);
    dh_sin(2:end) = diff(h_sin);
    K_alpha(n) = sum(dh_sin.^2)/(pi*dx^2);

    h_cos = cos(n*x_EW);
    dh_cos = h_cos;
    dh_cos(1) = h_cos(1)-h_cos(end);
    dh_cos(2:end) = diff(h_cos);
    K_beta(n)  = sum(dh_cos.^2)/(pi*dx^2);    
end
ModelSettings.M2ModeWeighting=[K_alpha,K_beta];

% Parameters about the model
% Assuming the data from KMC is in a domain [-L,L]
ModelSettings.MV      = 'T';
ModelSettings.FitTo   = 'meanR2';
ModelSettings.FittingTimeRange = 1000;
% ModelSettings.W       = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50];
ModelSettings.T       = [300,400,450,500,550,600,620,630,640,650,670,700];
ModelSettings.nu      = [2.6068e-003,2.5277e-003,2.4379e-003,1.0669e-003,0.4869e-003,0.2657e-003,9.6627e-005,8.0198e-006,0.8427e-003,6.6350e-003,1.855e-002,4.8547e-002];
ModelSettings.sigma2  = [0.6307,0.6179,0.6334,0.4737,0.3096,0.1757,0.0291,0.0028,0.0131,0.0395,0.0638,0.0971];
ModelSettings.Rh      = [0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2];
ModelSettings.K       = [0.6440 0.6439 0.6520 0.7676 0.9138 0.9823 0.9876 0.9902 0.9909 0.9916 0.9930 0.9939];
ModelSettings.Tau     = [3.9399 3.9226 4.0073 4.6479 3.5363 2.9746 2.7827 2.6868 2.6295 2.5721 2.4575 2.3667];

ModelSettings.h = 0;
ModelSettings.rho = 0;
ModelSettings.meanAlpha2 = zeros(ModelSettings.mode,1);
ModelSettings.meanBeta2  = zeros(ModelSettings.mode,1);
ModelSettings.varAlpha2  = zeros(ModelSettings.mode,1);
ModelSettings.varBeta2   = zeros(ModelSettings.mode,1);
ModelSettings.meanR2     = model_CalMeanR2(ModelSettings);
ModelSettings.meanM2     = model_CalMeanM2(ModelSettings);

save model_T_001.mat ModelSettings

%%
T = [300,400,500,600,700]';
t=0:5:100;

H   = zeros(length(T),length(t));
SOR = zeros(length(T),length(t));
R2 = zeros(length(T),length(t));
M2 = zeros(length(T),length(t));
for i = 1:length(T)
    model = ModelSettings;
    for j = 1:length(t)
        model = Model_T_update(model,T(i),5);
        SOR(i,j) = model.rho;
        M2(i,j) = model.meanM2;
        R2(i,j) = model.meanR2;
        H(i,j)  = model.h;
    end
end
%%
figure(1);plot(t',H');xlabel('time (s)');ylabel('H');legend(num2str(T));
figure(2);plot(t',SOR');xlabel('time (s)');ylabel('SOR');legend(num2str(T));
figure(3);plot(t',R2');xlabel('time (s)');ylabel('R^2');legend(num2str(T));
figure(4);plot(t',M2');xlabel('time (s)');ylabel('M^2');legend(num2str(T));

%%
